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Abstract: We analytically compute the spectrum of the 
Hessian of the Hamiltonian for a system of N particles in- 
teracting via a purely repulsive potential in one dimension. 
Our approach is valid in the low density regime, where we 
compute the exact spectrum also in the localized sector. We 
finally perform a numerical analysis of the localization prop- 
erties of the eigenfunctions. 

Great efforts have recently been done to find a com- 
prehensive microscopic theory able to describe the be- 
haviour of super-cooled liquids Among the others, 
a key point is to understand the mechanism of the glass 
transition, and to find a description suitable both for the 
super-cooled liquid and for the amorphous solid which 
forms at low temperatures. 

The general frame where a great number of recent anal- 
ysis have been performed is the Instantaneous Normal 
Modes (INM) approach |2). The main idea of this ap- 
proach is that liquids are 'solid-like' at short times t < r, 
where the typical diffusion time r increases strongly when 
lowering the temperature. Liquids' dynamics would cor- 
respond in this picture to vibrations about some equilib- 
rium positions with periodic jumps into new local minima 
In order to describe in a quantitative way this be- 
haviour, it is important to study the properties of the 
Hessian matrix of the Hamiltonian, averaged over the 
equilibrium Boltzmann distribution. The crucial quan- 
tity is the typical spectrum of the Hessian, whose eigen- 
vectors are the so-called Instantaneous Normal Modes. 
The determination of the spectrum therefore represents 
a vital task for any theoretical study of liquids. 

The aim of this letter is to outline an analytic approach 
for the computation of the INM spectrum. Some im- 
portant steps in this direction have been recently done 
in |Q and J5j, where a remarkably good agreement with 
numerical simulations on liquids has been achieved and 
where more realistic models than the one we study here 
were considered. However, these former analytic compu- 
tations of the INM spectrum assumed a Gaussian proba- 
bility distribution for the auxiliary degrees of freedom of 



the liquid B (this technical point will be clarified later) 
and we think that a Gaussian approach is likely to be 
too simple for a study of non-trivial spectral properties, 
such as localization. This is a crucial difference with the 
computation we present here, which, going beyond the 
Gaussian approximation, allows us to better investigate 
the spectrum in the localized sector. 

Let us consider a system of N interacting particles with 
Hamiltonian H = J2^ >1 V(ru), where V(r) is a two-body 
potential. The Hessian matrix A is defined by, — 
didYH, with fi, v = 1, . . . , d, being d the dimension of 
the space. In the present letter we shall consider the 
one dimensional case d = 1, in order to keep the algebra 
as simple as possible. The form of A is Aki — —Jki + 
Ski Jki, where Jki = V"(r k i). The diagonal term of 
A is a consequence of the translational invariance of the 
system, which requires Am — 0. 

As a first step in our analysis, we must find out what 
is the probability distribution of the matrix J induced by 
the equilibrium probability over the positions of the par- 
ticles. We introduce here our first approximation, assum- 
ing that the probability distribution P[3] is factorized 
into the single probabilities of the particles pairs. In this 
way the elements of J (but not of A) are independently 
distributed, i.e. P[3] = Ylk>lP(^kl)- This approximation 
works well at low densities, where the three-particles cor- 
relations are negligible. 

Once assumed this factorized form for P[3], we can 
express the pair-probability p( J) as, 



P{J) 



N 



drgW(r) 5(J-V"(r)) 



(1) 
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where (r) is the two-particles correlation function and 
p is the density. In the following we shall consider the 
soft-spheres potential V(r) = l/r m . In this case den- 
sity and temperature scale homogeneously, so that we 
can simply fix the density p — 1 and consider the high 
temperature regime as equivalent to the low density one. 
This does not hold for a non-homogeneous potential and 
what follows must be interpreted as a low density calcu- 
lation. 

In the definition of p( J) we must insert the explicit 
form of (7^ 2 '(r), which can be obtained by means of some 
approximation schemes in liquid theory or by numerical 
simulations. In the present context we want to show what 
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are the results of our method in the simplest analytical 
way, so that we stop at the first order of the virial ex- 
pansion and assume g^ 2 '(r) = exp(— (3V(r)). Thus, from 
eq.(Q) we have, 



P(J) 



1 

N J 1+c 



N 



q{J) 



(2) 



with f3 s — (3 [m(m + 1)] m + 2 , b = m/(m + 2) and 
c = l/(m + 2). For realistic values of m (typically 
rn = 12) the parameter b is very close to one. Therefore, 
we will directly set b = 1 in p( J) in order to simplify our 
calculation. We will show in the discussion of our results 
that the actual spectrum is very weakly dependent on 
this approximation. As it stands the distribution p(J) is 
not normalizable, but we can regularize it in the follow- 
ing way. Let us put an IR cut-off f , by setting V(r) =0 
for r > f, and let r\ = V(f). We obtain in this way a 
regularized form of the pair probability: 



df q(J') 



(3) 



where q(J) is defined in equation 
we shall use the notation, 



In the following 



(f(J)) 



dJ q(J) {f (J) - / (0)] 



(4) 



Note an important point: the distribution Pri(J) is di- 
luted, since the probability of finding an element of the 
matrix J larger than r\ is of order l/N. 

In order to compute the spectrum of A we introduce 
the resolvent, Gfcz(A|J) = [A — A + ie] kl . In the follow- 
ing we shall include the small imaginary term ie in A. 
The spectrum (or density of states) D(X) is then given 
by, D(X) = lim^olm Tr G(A)/ttA, where G(A) is the 
average over J of the resolvent. Following the recursive 
method of [Q, we write now a self-consistent equation 
for G. Given a system of N particles, we add an extra 
particle, with label 0, and we write the expression of the 
(0, 0) component of G in the (N + l)-particles system, 
isolating the contribution of the 0-th particle: 



^ +1) (AU) = \ J 



dtp dipi . . . dtpN (f> 
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where Z is the effective partition function and where (j> = 
ipQ. The action S is given by, 

. JV . N 

s - W - \ E J °*(<t> - to? + 5 E MG (N) \- k ^ ■ 
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k,l=l 



In the same spirit as the cavity method || , we have writ- 
ten a self-consistent equation which relates the element 



G 



(N+l) 



trix Gjy of the system with TV particles. We note that 



of the system with (TV + 1) particles, to the ma- 



G^P is the resolvent of the system in absence of the 
particle and therefore it does not depend on the new vec- 
tor Jofc. Besides, it is shown in J?J that only the diagonal 
elements contribute to the last sum in the action S. 

Since we have to perform an average over J, it is useful 
to introduce replicas, by writing Z~ x — lirn n _ > o ■ 
After replication we have: 



iG, 
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(MJok,G^>) 



lim — 
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d<H^) 2 r^|j 0fe ,4f), 



where 4> = (</>i . . . 4> n ) and 

n($\J 0k ,Gf k ) ) = e i ^^ x 

N 



fc=i 



The crucial quantity is now the average probability distri- 
bution fl((f>), obtained by averaging expression (||) over J. 
Indeed, from f2 we can reconstruct the average resolvent 
G by integration over <p and therefore obtain the spec- 
trum. Thus, unlike what has been done in |7j, we shift 
our attention on the pursue of a self-consistent equation 
for Q, and not for the resolvent G itself. Let us note 
that it is this very distribution f2(</>) that, as said in the 
introduction, has been assumed to be Gaussian in the cal- 
culations of Q and H (the distribution s(X) in Section 
III.B of H is exactly the same object as our f2(0)). 

The dependence on the matrix J of Q is divided in two 
independent parts: the vector J$u and the matrix G^). 
We can therefore compute separately the two averages 
and, after some algebra, we have: 



exp 



We can now obtain a more tractable integral equation by 
defining the following function g: 

g {4>) = ( ( e -5 J 5::w«-*-) ! ) . 



Note that g measures the deviation from Gaussianity of 
f2, so that obtaining a non-quadratic form of g means 
going beyond the Gaussian approximation of Q and ||. 
The integral self-consistent equation for g is, 



g{4>) 



We assume now that g depends only on the modulus x 
of the replica-vector <f>, so that we can perform the an- 
gular part of the integral, then average over the disorder 
distribution p v (J) and finally let 77 — > 0, to get: 



r(2 - c) 

c(c-l) 



IX 



x \ dyK{x,y) jW+ato) 



(7) 
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dt 

o Vl-t 



r(3=c)r(3±c) 4 



These two equations are the main result of this letter. 

for 



It is possible to prove analytically that g(x) ~ — x 2c 
x 3> 1, thus proving that f2 is definitely noi a Gaussian 
distribution. Notably, we have been able to numerically 
solve the equation for oix) without any further approx- 
imation. Indeed, eq.(^) has the form of a fixed-point 
equation, so that it is tempting trying to solve it numeri- 
cally by iteration. This is what we have done, discretizing 
the function g and the kernel if on a lattice. We have 
found that the convergence is rather fast and very weakly 
dependent on the small imaginary part e of A. Indeed, 
by setting directly e = 0, the results are very satisfactory 
and we are able to obtain the solution up to arbitrary 
precision. Once obtained g for a given value of A we have 
computed the spectrum D(X), using the formula, 



D(X) = Re i - 

TT 



dx 



Jfx +g(x) 



which follows from the definitions of D, G and g. The 
results are shown in Fig.l, where we have plotted the 
INM spectrum D as a function of A, for m = 12. 
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FIG. 1. The INM spectrum D as a function of A for differ- 
ent values of the scaled temperature j3 s ; m — 12 and e = 0. 
The plot is in log-linear scale. Inset: D(X) for (3 S = 1.5 in 
linear-linear scale. The spectrum vanishes at A = 0. 



The spectrum has positive support because d — 1 and 
it depends on the scaled inverse temperature (3 S in the 
expected way: for low temperatures (high f3 s ) the colli- 
sions among particles are weaker, so that the spectrum 
is peaked on lower value of the eigenvalues. On the other 
hand, the tail for large A is larger at higher temperature. 



not been able to express a as a function of the parame- 
ters (3 and m, even if it is numerically evident that a is 
a monotonically increasing function of j3 8 . 

A crucial task is now to check whether the result we 
have found is correct. To this aim we have done ex- 
tensive numerical simulations. Once drawn a matrix J 
with probability (^|), we build A and diagonalize it nu- 
merically. Since the spectrum has huge tails for large 
eigenvalues, it is convenient in order to compare simu- 
lations with analytic results to consider the probability 
distribution tt of /i = In A, that is 7r(/i) = D(e /1 )e'\ In 
Fig. 2 we plot 7r(/i) as obtained from the analytic form 
of -D(A), together with the one obtained from numerical 
simulations. The two curves are in excellent agreement 
confirming the validity of our result. Besides, we show in 
the inset of Fig. 2 the numerical spectrum obtained with 
the original value of b = m/(m + 2). The result justify 
the sensibility of the approximation b ~ 1. 
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FIG. 2. Numerical simulations vs. analytic solution. We 
plot here for 6=1 the probability distribution 7r(/i), with 
ju = In A; N = 600, r] = 10" 4 , f3 s = 1 and m = 12. Inset: 
on the same scale, analytic result for b = 1 vs. simulations 
performed with b = m/(m + 2). 



An equation similar to (Jfy have been derived in ||, 
with a different method and within a different con- 
text. Also in that case the distribution of the disor- 
der was diluted and translationally invariant, but the ex- 
plicit probability distribution p(J) was bimodal, p{ J) = 
6 (J) + j$[S(J — 1/p) — S(J)], where p was the connectiv- 
ity ||. Yet, even with this simple bimodal distribution 
the spectrum had not be worked out explicitly until now. 
Remarkably, we have been able to compute exactly the 
spectrum associated to this bimodal distribution, by nu- 
merically solving the corresponding self-consistent equa- 
tion ||. The spectrum is shown in Fig. 3 and can be 
compared with the approximated solution of fl0|| . Note 
that in the liquid spectrum there is no trace of the small 
tails oscillations present in the bimodal case |nj . 
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FIG. 3. Analytic solution. The spectrum in the case of the 
bimodal distribution, with p = 20 and e = 0.005. The small 
non-zero value of D for A — > is due to the non-zero value of 
e, since the limit e — + is very delicate in the bimodal case. 

An important issue is the analysis of the localization 
properties of the eigenfunctions. With respect to this 
an important quantity is the average inverse participa- 
tion ratio (IPR) ^(A), which provides information on 
the nature of the eigenfunctions and can be easily com- 
puted via numerical simulations. The IPR is defined as, 
Y(\ a ) — ( w a) ' where a = 1...N is the eigen- 

value index and w l a = [(Aq,|z)] 2 is the weight of site i in 
the eigenfunction |A Q ). 



Y — > 1 for A — > is smooth, due to the fact that the J 
distribution we consider is divergent at J = 0, so that 
it never happens that one single element in a row of J 
is of order one, with all the others equal to zero. This 
behaviour in the left tail cannot be present in the case of 
a bimodal distribution. On the other hand, the localized 
states of the right tail correspond to pairs of very strongly 
interacting particles and this naturally leads to an IPR 
equal to 1/2. A more detailed discussion of localization, 
following |ll[, will be presented elsewhere. 

In this letter we have outlined a general method to 
study analytically the INM spectrum of a liquid at low 
densities. We have succeeded in an exact computation in 
the simple d = 1 case, but our method is suitable to be 
extended to dimensions larger than one, where the Hes- 
sian is not positive defined and negative eigenvalues exist. 
The presence on negative modes in three-dimensional sys- 
tems is particularly relevant in connection with the glass 
transition 0]. Indeed, it has been argued in the context 
of the INM approach, that the Mode Coupling transition, 
marking the crossover from a non-Arrhenius behaviour of 
the viscosity to an Arrhenius one, occurs when the frac- 
tion of negative delocalized modes (the only ones related 
to collective particles diffusion) drops to zero || . In view 
of this, an approach as the one presented here, able to 
investigate also the properties of the localized modes, can 
prove extremely useful. We will address the extension of 
our method to higher dimensions in a future work. 

The work of AC and IG was supported by EPSRC 
Grant No. GR/K97783. 
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FIG. 4. The IPR Y as a function of the eigenvalue A at 
different values of N. n — 10~ 4 , /3 S = 1 and m — 12. The 
thick curve is the corresponding spectrum D(X). 



In Fig. 4 we plot Y as a function of A. It is clear that 
there are two localizations edges, separating a central 
region of extended eigenvalues from the tails where lo- 
calized states are present. For A — ► the IPR goes to 
one and this corresponds to a single particle nearly de- 
coupled from the rest of the system. Note that the limit 
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